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Abstract 

Crop model-specific biases are a key uncertainty affecting our understanding of climate change impacts to agricul- 
ture. There is increasing research focus on intermodel variation, but comparisons between mechanistic (MMs) and 
empirical models (EMs) are rare despite both being used widely in this field. We combined MMs and EMs to project 
future (2055) changes in the potential distribution (suitability) and productivity of maize and spring wheat in South 
Africa under 18 downscaled climate scenarios (9 models run under 2 emissions scenarios). EMs projected larger yield 
losses or smaller gains than MMs. The EMs' median-projected maize and wheat yield changes were —3.6% and 6.2%, 
respectively, compared to 6.5% and 15.2% for the MM. The EM projected a 10% reduction in the potential maize 
growing area, where the MM projected a 9% gain. Both models showed increases in the potential spring wheat pro- 
duction region (EM = 48%, MM = 20%), but these results were more equivocal because both models (particularly the 
EM) substantially overestimated the extent of current suitability. The substantial water-use efficiency gains simulated 
by the MMs under elevated CO 2 accounted for much of the EM— MM difference, but EMs may have more accurately 
represented crop temperature sensitivities. Our results align with earlier studies showing that EMs may show larger 
climate change losses than MMs. Crop forecasting efforts should expand to include EM— MM comparisons to provide 
a fuller picture of crop-climate response uncertainties. 
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Introduction 

Climate change impacts focusing on crops should 
account for several important sources of uncertainty 
when quantifying the range in potential future effects. 
The most obvious uncertainty is the extent to which 
regional precipitation and temperature patterns will 
change. Most studies address these issues using inputs 
from climate model ensembles that capture the range of 
variation in climate trajectories caused by climate 
model biases and possible emissions trajectories. 
Another important source of uncertainty is the struc- 
tural differences between crop models, which has 
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started to become an active research focus in the past 
few years. Crop model uncertainty is the major focus of 
the Agricultural Model Intercomparison and Improve- 
ment Project (AgMlP, Rosenzweig et al., 2013), and has 
been the focus of several studies (e.g., Palosuo et al., 
2011; Tubiello et al., 2007; Rosenzweig & Wilbanks, 
2010; Rdtter et al., 2011; Knox et al., 2012). 

However, despite this growing awareness of inter- 
model variation among agricultural impacts research- 
ers, little attention is given to the potentially larger 
differences between empirical and mechanistic (or pro- 
cess-based) modeling approaches. Multimodel compar- 
isons for crop forecasting compare only mechanistic 
models (MMs) (e.g., Asseng et al., 2013; Rickebusch 
et al., 2008; Rdtter et al., 2011; Challinor & Wheeler, 
2008; Palosuo et al., 2011), yet both empirical and MMs 
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have strengths and weaknesses that should be taken 
into account. 

Ecological impacts researchers have explicitly com- 
pared empirical models (EMs) and MMs, providing 
insight into how these two classes differ (Dormann, 
2007; Morin & Thuiller, 2009; Cheaib et al., 2012; 
Dormann et al, 2012; Hijmans & Graham, 2006). These 
studies found that EMs tend to project larger plant 
range shifts than MMs (Morin & Thuiller, 2009; Cheaib 
et al., 2012), which may be due to EMs' tendency to 
underestimate persistence under novel climates, or 
poor ability to simulate plant responses to elevated 
CO 2 (Morin & Thuiller, 2009; Cheaib et ah, 2012). On 
the other hand, MMs typically do not simulate 
potentially limiting biotic interactions (e.g., interspecific 
competition), and thus may understate negative 
impacts if warming exacerbates such constraints 
(Morin & Thuiller, 2009; Tubiello et al, 2007). EMs 
implicitly incorporate these factors because training 
data reflect the species' current distribution, which is 
partly shaped by biotic constraints (Morin & Thuiller, 
2009; Cheaib et al, 2012). 

Empirical models thus represent the species' perfor- 
mance within its 'realized' environment, which includes 
all biotic and abiotic limits, whereas MMs characterize 
potential performance in response primarily to abiotic 
factors (Cheaib et al, 2012). For crop species, EMs simu- 
late realized yields in response to prevailing farm 
conditions (current management, pests, diseases, soil, 
water, climate), while MMs simulate how potential 
yields (given existing abiotic conditions) may vary in 
response to management practices. Given these struc- 
tural differences, crop EMs and MMs should show 
climate change responses similar to those found in 
ecological studies, with EMs showing larger yield (and 
potential growing region) losses than MMs. Direct com- 
parisons of crops' climate responses simulated by inde- 
pendently developed EMs and MMs are largely absent 
(but see Lobell & Burke, 2010; for a comparison of EMs 
fit to MM-generated data). To date, only a single study 
of maize compared crop suitability and yield predicted 
by different modeling approaches, but this work did 
not project future climate impacts (Estes et al, 2013). 

Systematic differences between crop EMs and MMs 
are apparent when comparing results from different 
agricultural impacts studies. For example, results for 
South Africa from several regional-scale impacts stud- 
ies indicate that maize yields will decline 10-30% by 
2050 (Parry et al., 1999; Jones & Thornton, 2003; Schlen- 
ker & Lobell, 2010; Knox et al, 2012), wheat yields will 
be 16% less by 2030 (Lobell et al, 2008), and overall 
cereal productivity will fall 0-50% by 2100 (Parry et al, 
2004; Fischer et al, 2005). Among these studies, the two 
using EMs (Lobell et al, 2008; Schlenker & Lobell, 2010) 


projected yield losses that were twice as large (28-30%) 
as the two most comparable MM-based results (10-19%, 
Jones et al, 2003; Parry et al, 1999; summarized by 
Schlenker & Lobell, 2010) . Although these differences 
may be partly due to the low fertilizer inputs used in 
the MM studies (which would have reduced projected 
climate sensitivity; Schlenker & Lobell, 2010), given the 
food security implications of these findings, and the 
aforementioned evidence regarding EM-MM discrep- 
ancies from the ecological literature, studies directly 
comparing agricultural impacts projections from these 
two model classes are needed. 

In this study, we compared EM and MM projections 
of climate change impacts to maize and spring wheat, 
two important cereals that use different photosynthetic 
pathways (C4 for maize; C3 for wheat). We focused 
on South Africa, a regionally and globally important 
agricultural power that is the world's 9th largest maize 
producer and sub-Saharan Africa's 2nd largest wheat 
grower (FAO, 2012). Although it is semiarid, most 
(85%) of South African maize and wheat is rainfed 
(Hardy et al, 2011). Some climate scenarios show that 
South Africa will become both drier and hotter (IPCC, 
2007), which could cause substantial production losses. 
We undertook a spatially explicit investigation of how 
the potential growing areas and productivity of rainfed 
maize and wheat might be impacted by mid-21st cen- 
tury climate change. Our goal was to examine how crop 
impacts projections from EMs and MMs differ, while 
providing further insight into how climate change may 
impact regional food security. 

Data and methods 

Background and model overview 

South African maize is grown primarily in the summer rain- 
fall (ca. 80% of precipitation between October and April) 
region in South Africa's northeastern half (Fig. 1), while the 
bulk of wheat (ca. 50%) is produced in the year-round and 
winter rainfall regions (ca. 55-80% of rainfall between April 
and September) along the southern and southwestern coasts 
(Hardy et al, 2011). Maize plantings and production average 
31 000 km^ and 10 000 kt, respectively, while wheat extent 
and production is 8200 km^ and 2300 kt, respectively (Anony- 
mous, 2009). Less than 10% of either crop's production area is 
irrigated (Bradley et al, 2012). 

We simulated future changes in maize and spring wheat 
yields, as well as shifts in these crops' potential growing 
regions (hereafter 'suitability'). We used a widely employed 
mechanistic crop growth simulator (MM) and a simpler EM 
run with climatological variables. The MM we selected was 
version 4.5.0.047 of the Decision Support System for Agrotech- 
nology Transfer Cropping System Module (DSSAT-CSM, 
hereafter 'DSSAT' Jones et al., 2003; Hoogenboom et al, 2012), 
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Fig. 1 South African maize and wheat-growing regions. Gray 
points denote aerial observations (collected between 2006 and 
2009) of crop fields growing maize, while black points are fields 
growing spring wheat (2006-2008). 


which incorporates the well-established Crop Environment 
Resource Synthesis (CERES) maize and wheat models (Ritchie, 
1998). DSSAT models simulate growth and reproduction in 
relation to environmental conditions and management prac- 
tices (e.g., fertilizer applications, planting densities, cultivar 
characteristics) at a daily time step (Jones et al., 2003). CERES 
models simulate potential daily growth (PG) rates using a 
radiation-use efficiency (RUE) approach: 


PG = 


RUE X PAR 
n 


(1 _ e~'‘^^)C02 


( 1 ) 


where PAR is photosynthetically active radiation 
(MJ m“^ day^'), n is plant density per m^, fc is a light extinc- 
tion coefficient, LAI is the leaf area index, and CO 2 scales 
growth according to atmospheric CO 2 concentrations (Hoo- 
genboom et al., 2012; Ritchie, 1998). The growth rate is further 
modified by factors (calculated in submodules) related to 
nitrogen and water stress and temperature and soil fertility, 
and dry matter is partitioned according to crop growth stage 
(determined by growing degree days). The CO 2 factor, which 
is the ratio of photosynthetic rates under elevated CO 2 relative 
to reference CO 2 levels, was initially derived from chamber 
experiments, but later updated to reflect Free Air Carbon 
Enrichment (FACE) experiment results (Rosenzweig & Igle- 
sias, 1998; Backlund et ah, 2008; Boote et ah, 2011; Hoogen- 
boom et ah, 2012). DSSAT's evapotranspiration module 
increases water-use efficiency under higher CO 2 (via stomatal 
closure) using a factor derived from the ratios of elevated CO 2 
to reference transpiration rates (Rosenzweig Iglesias, 1998). 

We used Generalized Additive Models (GAMs; Yee (& 
Mitchell, 1991) to develop EMs for maize and wheat. The 
GAM is a variant of the generalized linear model that can 
accommodate complex nonparametric relationships between 
each predictor and the dependent variable, and takes the 
form: 


g(y) = Po + Si (Xl,) + S2(x2i) + . . . + (2) 

where g is a link function allowing nonlinear relationships 
between the dependent variable y and predictors Xi_j, Pq is 
the intercept parameter, and Si_y are optional smoothing func- 
tions that can be replaced by linear terms (in which case s is 
replaced by p Yee & Mitchell, 1991). GAMs are widely used in 
ecology, can be fit to both binary and continuous response 
variables, and perform well in estimating current South Afri- 
can maize productivity and distribution (Estes et ah, 2013). 

In this study, we examined the relative changes in produc- 
tivity and suitability between the periods 1979-1999 ('the base- 
line') and 2046-2065 ('the future'), as projected by the two 
different models developed for each crop, based on input 
from an ensemble of downscaled climate scenarios. 

Datasets 

Crop productivity and distribution data. To develop and val- 
idate our models, we used two remote sensing datasets. The 
first was airborne crop census data from 2006 to 2009, together 
with accompanying crop field boundary maps (Eourie, 2009; 
SiQ, 2007). The second was 16 day normahzed difference veg- 
etation index (NDVI) imagery collected by the Moderate Reso- 
lution Imaging Spectrometer (MODIS; source: http:/ /Ipdaac. 
usgs.gov) for 2006-2009. These two datasets, respectively, pro- 
vided distribution points (14 736 for maize; 1355 for wheat) 
and proxy yield variables for the two crops, the development 
of which is described by Estes et ah (2013) for maize. We 
summarize these methods here (see Appendix SI) and 
extended them to spring wheat. We used the crop census 
points to select NDVI time series (for the year the census 
observation was made) that were primarily composed of 
either maize or wheat reflectances, and integrated the NDVI 
values over the growing season to create a remotely sensed 
estimate of each crop's yield ('NDVIyieid'). After averaging 
points to a 20 km resolution (to minimize site-specific and in- 
terannual variability), there were 436 NDVlyi^ij values for 
maize and 122 for wheat. 

We developed the NDVIyi^u dataset because reported 
yields were only available for South Africa's provinces, which 
have a mean area of 135 646 km^. Although the provincial 
data could have been used in a panel modeling approach to 
estimate crops' climate responses (by incorporating both the 
interannual and the interprovincial variation in yields; Lobell 
& Burke, 2010), we wanted to model yields at a scale that was 
closer to the resolution of our input soil and baseline climate 
datasets (<135 km^). We thus used NDVIyieid, which (when 
aggregated to the provincial scale) explains 59-67% of 
variance in reported yields. Other studies also confirm that 
this measure is an effective maize yield proxy (see Estes et ah, 
2013 for citations and Appendix SI). For wheat, there 
were insufficient observations to robustly assess the proxy 
against provincial reported yields (Appendix SI), but NDVl- 
derived wheat yield estimates have been effective for other 
regions (e.g., = 0.64 in Wang et ah, 2005; and see Ren et ah, 

2008). 

The 20 km aggregates of NDVIyieid were used for training 
and testing the GAM maize and wheat models, and to assess 
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the accuracy of DSSAT baseline yields. The crop census data 
were used to create the crop suitability maps (see below). 

Soil data. Decision Support System for Agrotechnology 
Transfer requires data describing soil drainage rate, horizon 
depth, wilting point (WP), field capacity (FC), saturated mois- 
ture content, bulk density, and organic carbon content. Estes 
et al. (2013) mapped these values by applying pedotransfer 
functions to soil texture and composition data from South 
Africa's Land Type database (SIRI, 1987), which we used here 
to run DSSAT maize and wheat. We used a subset of these to 
fit GAM models [FC, WP, plant available water (FC minus 
WP), soil depth, and organic carbon; see Appendix SI]. 

Climate data. For the baseline, we obtained daily values of 
rainfall, minimum and maximum temperature (Tmin/ Tmax)/ 
and solar radiation from the South African Quinary Catch- 
ment database (Schulze & Horan, 2010), which has a 50 year, 
spatially and temporally rectified weather record linked to 
5838 river basins (mean area 135 km^) covering South Africa, 
Lesotho, and Swaziland. 

For the future period, we used an ensemble of 9 Coupled 
General Circulation Models (CGCMs; see Appendix S2 for 
list) run under both a relatively high (A2) and low (Bl) emis- 
sions scenario (Nakicenovic & Swart, 2000). This allowed us 
to account for uncertainties related to the future emissions 
trajectory as well as the simulated atmospheric response. 
CGCMs provide coarse resolution (200-300 km) climate 
change estimates, but crop impacts studies typically require 
higher resolutions. We downscaled the CGCMs to the Qui- 
nary Catchment resolution using an empirical method based 
on self-organizing maps (Hewitson & Crane, 2006). This tech- 
nique (summarized in Appendix SI) provides realistic precip- 
itation values that are more consistent with observed patterns 
in South Africa than those generated by the 'parent' CGCMs 
(Hewitson & Crane, 2006), and was shown to provide consis- 
tent results between CGCMs across Africa in the IPCC 4th 
assessment report (Christensen et al., 2007). The method has 
been used extensively to downscale climate in southern 
Africa (Crespo et al., 2010; Patt et al., 2010; Tadross et al., 
2009). 

While the resulting rainfall simulations were more realistic 
and consistent (Hewitson & Crane, 2006), the downscaled cli- 
mate data were still subject to CGCM-specific biases. To 
remove the remaining discrepancies between the observed 
basehne weather records and the downscaled CGCM base- 
lines, we adjusted the distribution parameters of the daily 
records. We first calculated the monthly climatological anoma- 
lies between each downscaled CGCM baseline-future pair, 
and used these to adjust the distribution parameters for Tmin/ 
Tmax/ arid daily rainfall, as calculated for each month in the 
20 year baseline record (i.e., the climatologies of the observed 
basehne). We altered the mean and SD of Tmm and Tmax/ and 
the mean and the shape and scale parameters of the gamma 
distribution fit to observed rain days. We also altered the 
number of rain days to match the change in rainfall frequency 
reflected in each CGCM basehne-future pair. The shape of the 
cumulative distribution function for each variable for each 


month was then iteratively adjusted until its summary statis- 
tics matched the adjusted statistics, within a specified error 
tolerance (0.3 °C for mean Tmax arid Tmm and 0.05 mm for 
mean daily rainfall, 0.1 for SDs and 0.075 for gamma shape 
and scale). 

The procedure transformed the baseline weather data into 
bias-corrected, future daily values for Tmm, Tmax, and rainfall 
for each of the 18 downscaled CGCMs. Solar radiation was 
unadjusted except where rain days were adjusted, in which 
case it was either reduced by 10% (basehne dry days con- 
verted to future rain days) or increased by 10% (baseline rain 
days converted to future dry days). We evaluated the proce- 
dure's accuracy based on the root mean squared error 
between the summary statistics derived from the (i) monthly 
climatologies of the bias-corrected future daily records; 
and (ii) the observed baseline climatologies that were initially 
adjusted by the CGCM control-future anomalies (see Appen- 
dix SI). 

The bias-corrected daily records were used to run DSSAT. 
For the GAMs, we calculated the following climatological 
variables from the bias-corrected daily data: average cumula- 
tive growing season precipitation, the mean maximum Tmax 
and the mean minimum Tmm of the coldest and hottest 
months in the growing season (where seasons 243 are Octo- 
ber-April for maize and May-September for wheat Bradley 
et al., 2012; Estes et al., 2013). The median projected changes in 
growing season precipitation (both summer and winter) and 
maximum Tmax are illustrated in Fig. 2. The downscaled pre- 
cipitation results contrast with the IPCC AR4 ensemble mean, 
which projected reduced daily mean precipitation for the 
whole country by 2080 (IPCC, 2007). 

Model development and validation 

Mechanistic models. The DSSAT model simulates crop 
responses to both environmental (described above) and 
management variables. The key management variables 
relate to planting (row spacing, plant density, sowing date), 
fertilization, and cultivar selection (Jones et al., 2003). We 
used industry recommendations and national cultivar trial 
data to define these parameters for both crops. Estes et al. 
(2013) describes the management parameter development 
for maize, and the sources and derivation steps for wheat 
parameters are detailed in Appendix S2. Planting density 
for both crops was determined using simple linear relation- 
ships found between mean annual rainfall (MAR) and data 
from the aforementioned source material. Row spacing was 
fixed at 25 cm for wheat and varied with rainfall for maize 
(between 0.9-2.1 m). Planting occurred automatically for 
both crops within specified date ranges when plant avail- 
able moisture exceeded 70%. The date range was May 5th 
to June 25th for wheat. For maize, we determined an aver- 
age planting date linked to MAR (from cultivar trials), and 
automatic planting occurred up to 2 weeks prior to this 
date and as late as January 15th. For both crops, 32 kg ha^' 
of Nitrogen fertilizer was applied at planting. For maize, 
we also ran the model with 59 kg ha^' N to test how 
higher rates impact model performance, given that average 
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fertilization practices are potentially higher than the selected 
level (see Appendix S2 for N input details). 

The DSSAT model uses coefficients to represent the envi- 
ronmental responses of different crop cultivars (Jones et al., 
2003). Ideally, these coefficients should undergo site-specific 
calibration, but in practice, the required data are rarely avail- 
able. We applied previously used maize coefficients (Estes 
et ah, 2013) representing a generic short-medium season 
cultivar. For wheat, we selected coefficients for an Austrahan 
cultivar (Appendix S2) grown in climates similar to that of 
South Africa's spring wheat-growing region. 

To run DSSAT spatially, it was necessary to define fields 
representing unique climate-soil combinations. Merging the 
soil and climate units (defined by the Quinary Catchment 
boundaries) created 107,140 'fields' with a mean area of 
1 1 km^. We ran simulations for both crops for each year in the 
1979-1999 period, averaged yields over this period, and aggre- 
gated results into 923 m resolution grids (to compare with the 
GAMs). We evaluated the accuracy of DSSAT's simulated 
yields based on the strength of their correlations with 
NDVIyieM- 

We constructed binary suitability maps for the two crops 
using thresholds of simulated yield and its coefficient of 
variation (CV, the SD of yields simulated over the 20 year 
period divided by the mean annual yield and converted %). 
Areas with mean yields above the yield threshold and yield 
CVs below the CV threshold were classed as suitable, while 
areas with one or both of these values on the other side of 
its threshold were unsuitable. Thresholds were found using 
the crop census observations and an equal number of 
randomly selected, pseudoabsence points placed outside 
the relevant crop growing regions to extract simulated yield 


and yield CV values. Thresholds for yield and yield CV 
were selected based on their ability to maximize the number 
of true positives (census points) and true negatives (pseud- 
oabsences). The accuracy of the resulting binary suitability 
maps was assessed against observed suitability surfaces 
created by creating kernel density estimates of mapped crop 
fields within the extent of each crop's growing area as 
defined by the aerial census data [see Estes et al. (2013) and 
Appendix S2 for further details on this and the calculation of 
thresholds]. 

Empirical models. We adapted the approach developed by 
Estes et al. (2013) to fit and validate GAMs for the two crops. 
We used NDVIyi^id as the response variable (n = 436 for 
maize, 122 for wheat), and the predictors were mean total 
growing season precipitation, and the average minimum and 
maximum temperatures of the coldest and hottest month in 
the growing season (October- April for maize, May-September 
for wheat Bradley et al., 2012; Estes et al., 2013). To allow for 
more realistic drought sensitivity, we also included the soil 
variables described above. In fitting GAMs, we applied 
smoothing functions to predictors exhibiting nonlinear rela- 
tionships with the response variable, provided that the 
smoothing function significantly improved the model fit 
(Wood, 2001) and had a biologically meaningful interpreta- 
tion. We discarded variables that were redundant or did not 
explain a significant proportion of variance in the response 
variable. 

We used 10-fold cross-validation to assess the accuracy 
of GAM-predicted yields and robustness of the selected 
model structures. To create suitability maps from the 
GAM-predicted maize and wheat yields, we used the same 
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Fig. 2 Baseline (1979-1999) climatologies and their median projected anomalies derived from the bias-corrected results of the 18 
downscaled climate scenarios. The left column shows the mean precipitation falling in the summer (October-April, top) and winter 
(May-September, bottom) halves of the year, and the second column displays median projected rainfall changes. The third column 
illustrates each season's mean maximum Tn,ax, and the fourth column shows the projected changes in this variable. 
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thresholding procedure described for DSSAT minus the 
yield CV threshold, as the GAMs did not produce annual 
yield estimates. 

Crop projections. Each crop model was run with the 18 
bias-corrected climate scenarios. DSSAT models can directly 
simulate the effect of differing CO 2 concentrations on crop 
growth and reproduction, whereas the GAMs cannot. For 
DSSAT, we set future CO 2 levels to 548 or 493 ppm, the 
respective means of the SRES A2 and B1 emissions scenarios 
for 2046-2065. To examine DSSAT's sensitivity to elevated 
CO 2 , we also ran each set of future scenarios with DSSAT's 
CO 2 inputs set to the baseline level. To test how much 
DSSAT's automatic (i.e., adaptive) planting methods affected 
crop yield projections, we ran each crop model using the cli- 
mate scenario that projected the median future rainfall for 
that crop's growing season, while keeping the planting dates 
fixed to those used in the basehne simulation. To assess how 
higher fertilizer application rates impact maize yields under 
future climates, we used the same climate scenario selected 
for the planting date analysis and ran DSSAT with 
59 kg ha^^ of N. 

We applied a post-hoc adjustment to GAM yield projec- 
tions to simulate increased photosynthetic and reduced tran- 
spiration rates under elevated GO 2 , following Pongratz et al. 
(2012). We increased GAM wheat yields by 13%, which we 
estimated using PAGE experiment results reported by Ains- 
worth et al. (2008) [we found the mean yield increase 
(14.25%) and the mean CO 2 increase (184 ppm) and rescaled 
this to the mean CO 2 increase for our study (173 ppm)]. For 
maize, FACE experiments suggest that it remains unclear 
how much elevated CO 2 enhances the photosynthetic rate of 
C4 crops, but reduced transpiration rates under higher CO 2 
may increase yields by lowering water stress (Leakey, 2009; 
Kimball, 2011). To account for this effect, we increased GAM- 
projected maize yields by 5% [proportional to Pongratz 
et al's. (2012) multiplier]. 

To evaluate changes in crop suitability, we created binary 
suitability maps from each simulated future yield (and yield 
CV, for DSSAT) surface. We found the fraction of suitability 
maps agreeing regarding future suitability (n/18), and used 
the median fraction (0.5) to assess the projected change in suit- 
ability. To assess climate impacts to productivity, we found 
the median future yield for each crop and calculated the per- 
cent change relative to baseline yield. We compared projected 
yield changes between models within the area shared by both 
models' baseline and future suitabilities. 

Results 

Baseline simulations 

DSSAT's simulated baseline maize yield explained 
37% of spatial variance in NDVIyjgid (Estes et al., 2013; 
Appendix S3). DSSAT maize yield aggregated to 
the provincial scale was strongly correlated with, but 
systematically underestimated, both mean reported 
provincial yields (R^ = 0.66) and aggregated NDVIyieid 


= 0.72; Appendix S3). Running the model with 
higher N (59 kg ha^') did not appreciably alter the 
relationship between simulated yields and NDVIyigid 
(R^ = 0.37), and resulted in a slightly poorer correla- 
tion with provincial yields (R^ = 0.58), although the 
bias toward yield underestimation was slightly 
reduced (Appendix S3). DSSAT's simulated baseline 
wheat yield explained 40% of spatial variance in 
NDVIyieid (Appendix S3). There were insufficient 
observations to assess provincial-scale correlations for 
DSSAT wheat. 

The best fitting GAM maize model included non- 
parametric relationships with mean summer precipita- 
tion and WP and linear relationships with mean 
maximum and soil depth. The selected wheat 

GAM model had nonparametric terms for average win- 
ter precipitation, mean minimum and WP, and a 
linear term for mean maximum Tmax- The maize GAM's 
10-fold cross-validation had average R^ values of 0.66 
and 0.64 for the training (the NDVIyieid values used to 
fit the model) and test partitions (the excluded 
NDVIyieid values), respectively, while the wheat GAM 
had mean R^ values of 0.73 and 0.57 for training and 
testing, respectively (Appendix S3). 

Both DSSAT and GAM's maize suitability showed 
similar total accuracy (true positives + false positives 
divided by total area) when compared to the observed 
maize suitability surface (87% and 86%; Appendix S3). 
The DSSAT suitability map was created using a yield 
threshold of 933 kg ha^' and a yield CV threshold of 
109%, while a yield threshold of 2957 kg ha^' was used 
for the GAM suitability map. 

The wheat suitability maps for the two models also 
showed comparable accuracy (95% for DSSAT and 93% 
for GAM; Appendix S3), but these statistics were 
inflated by the large size of the area unsuitable for 
spring wheat (115 000 km^) relative to the potential 
production region (5200 km^), and the fact that both 
models correctly classified most of the unsuitable area 
(negative predictive performance >0.98). However, both 
models substantially overestimated spring wheat suit- 
ability and thus had low positive predictive power (or 
precision; 0.46 for DSSAT, 0.33 for GAM). The GAM's 
false positive error extended up to the length of South 
Africa's east coast (Appendix S3). The suitability 
thresholds for DSSAT were 1150 kg ha^' for yield and 
70% for yield CV, and the GAM threshold was 
1433 kg ha^^ (Appendix S3). 

Projections 

Suitability. Both DSSAT and GAM projected that the 
core maize growing region will remain suitable in 2055 
(Fig. 3), based on the criterion of >50% model 
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agreement. However, the direction of suitability 
changes differed between models: DSSAT projected a 
9% gain in suitability, primarily along the southwestern 
boundary of the baseline growing region, whereas the 
GAM showed a 10% loss concentrated along the wes- 
tern and northern boundaries (see Appendix S4 for 
suitability changes under different agreement frac- 
tions). 

Both models showed an expansion in the total area 
suitable for wheat, with the GAM projecting expansion 
twice as large as DSSAT's (48% vs. 20%, Fig. 3). There 
was little loss projected by either model for its baseline 
region of suitability, although GAM showed small 
contractions in the extreme northwest. DSSAT showed 
suitability expanding into the interior from the north- 
west and southeast of the current growing region, 
whereas GAM showed gains into the interior along the 
entire length of its baseline suitability region (Fig. 3). 

To assess the sensitivity of DSSAT's suitability pro- 
jections to CO 2 , the simulations run with future climate 
scenarios, but under baseline, atmospheric CO 2 concen- 
trations were also converted to suitability maps. For 
maize, the median future suitability fraction showed a 
4% suitability loss and for wheat, a 3% gain. DSSAT 
therefore showed positive CO 2 effects of 13% and 17% 
for maize and wheat, respectively, in response to a 
projected mean CO 2 increase of ca. 170 ppm (the mean 
projected increase across the A2 and B1 scenarios rela- 
tive to the mean baseline concentration; see Appendix 
S4 for further details). 

Projected productivity. The maize productivity patterns 
simulated by the two models were substantially differ- 
ent (Fig. 4). DSSAT's baseline yield maps showed 
marked variation in response primarily to rainfall, with 
the highest yields (>4500 kg ha^') in the wetter, cooler 
highlands (Fig. 2) southwest and northeast of Lesotho 
(Fig. 4). Yields were typically >50% lower throughout 
the rest of the growing region. This spatial pattern in 
future projected yields was largely unchanged, but 
there were yield gains through much of the potential 
growing region. Gains were most pronounced (exceed- 
ing 200%) to the South of Lesotho, and in the high 
yield areas along the escarpment to the north of Leso- 
tho, and extending from there to the northwestern inte- 
rior (Fig. 4). On the east coast, yield was projected to 
rise by ~25%, likely due to less warming (<1.8 °C 
increase in Tmax) coupled with increased rainfall 
(>41 mm; Fig. 2). Yield losses of up to 25% were pro- 
jected in a belt extending from the northwestern 
boundary of Lesotho to the westernmost extreme of the 
maize growing region, with several small patches 
shows losses of 25-50%. Smaller yield losses (not 
exceeding 10%) were projected for the band lying 


along Lesotho's eastern border, and to the northwest of 
Swaziland (Fig. 4). 

In contrast, GAM's spatial yield patterns, both in 
terms of simulated potential and percent change, were 
more uniform than DSSAT's. GAM's highest potential 
areas corresponded to DSSAT's, and were of similar 
magnitude (>4500 kg ha^'), but the high yield areas 
were larger, and there was smaller range in yields and 
a more normal distribution of yield classes (Fig. 4). The 
GAM model projected uniform yield losses of up to 
10% throughout most of the growing region, reaching 
as high as 25% in the extreme west and north. Small 
gains (<10%) were projected for a few patches along the 
east coast and on the escarpment near Lesotho. 

Comparing the two models' projected yield changes 
within their jointly suitable areas, DSSAT's median 
projected maize yield change was +6.5% (interquartile 
range = —2.3 to 18.1%), while GAM's was —3.6% 
(interquartile range = —6.0 to —2.1%). 

The baseline wheat productivity patterns projected 
by DSSAT and GAM generally agree, with both models 
showing the highest yields (>2600 kg ha^') in the cur- 
rent core of the spring wheat production region in the 
southwestern Cape, and to a lesser extent in a narrow 
strip along the southern coast (Fig. 5), the two regions 
with the greatest winter rainfall (Fig. 2). However, the 
two models' projected yield change patterns diverged 
considerably. DSSAT showed an eastwards shift in 
productivity, with yields increasing by 25% or more 
along the southern coast, and ranging between +10% 
and —10% on the west coast. The GAM projected small 
losses along the western and southern coasts balanced 
by strong increases (>50%) into the interior of these 
regions. The east coast where GAM falsely showed 
baseline suitability was also projected to show some 
gains (Fig. 5). In the two models' jointly suitable areas, 
DSSAT's median projected wheat yield change was 
+15.2% (interquartile range = +6.5 to 19.7%), while 
GAM's was +6.2% (interquartile range = —2.6 to 39.6%). 

We assessed the CO 2 sensitivity of DSSAT's projected 
yield change estimates by comparing the median yield 
change from simulations conducted under projected 
CO concentrations with the median of simulations run 
under baseline CO 2 levels. For maize, there was a posi- 
tive CO 2 effect of 15.7%, and for wheat, the effect was 
+26.8% (see Appendix S4 for further sensitivity results). 
For planting date sensitivity tests, we reran DSSAT 
models with MRI A2 scenario for maize and CSIRO B1 
for wheat (the respective median rainfall scenarios) 
with planting dates fixed to those used during the base- 
lines simulations, and compared yields with the initial 
projections within the area of baseline suitability. Maize 
planting occurred 5.6 days earlier when planting was 
automatic, but yields were 14.8 kg ha^' lower on 
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Fig. 3 The change in areas suitable for maize and spring wheat production in 2055, as projected by both the Decision Support System 
for Agrotechnology Transfer (DSSAT) and Generalized Additive Model (GAM) models under the median agreement criterion (at least 
9/18 simulations agree regarding future suitability). Plots along the bottom row show the percent change in suitable area (given in 
km^ X 1000) relative to each model's simulated baseline suitability. 


average compared with fixed planting dates. For wheat, 
planting was 0.6 days earlier, and yields were 
4.7 kg ha^' lower. For the maize N sensitivity test, we 
reran DSSAT maize under the MRI A2 scenario with 
59 kg ha^' N, which showed a median yield change 
(relative to its baseline) of +19%, compared to +13.8% 
for the same scenario run with 32 kg ha^'. The spatial 
pattern of yield change under higher N was the same or 
more positive in 82% of the maize growing region, and 
in a further 13% the projections were at most 3.1% lower 
than the original treatment's results (Appendix S4). 

Discussion 

Empirical vs. mechanistic models 

Similar to a previous assessment (Estes et al., 2013), the 
GAMs' baseline yield and suitability estimates were 


more accurate than DSSAT' s, except in delineating the 
spring wheat suitability region. The better performance 
and the lower input data requirements suggest that 
simple EMs such as GAMs can be better choices than 
data intensive MMs for studying current crop-climate 
relationships over large spatial extents (Estes et al., 
2013). 

However, in terms of understanding crop responses 
to changing climates, the GAMs consistently projected 
larger yield and suitability losses than DSSAT. DSSAT 
models projected net productivity and suitability gains, 
driven by increasing yields in the east of both crops' 
current growing regions (Figs 4 and 5). GAMs showed 
either modest yield losses or smaller yield gains than 
DSSAT in these areas. The exception to this pattern 
was GAM'S large spring wheat suitability gains (dri- 
ven by limited model calibration data, see Caveats sec- 
tion below), caused by large projected yield increases 
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Fig. 4 Decision Support System for Agrotechnology Transfer- (top row) and Generalized Additive Model (GAM)-simulated (bottom 
row) maize yields for the baseline (left column) and future (middle column) periods, and the differences (right column) between the 
two expressed as a percentage. 


in the interior outside the baseline suitability region 
(Fig. 5). Along the southern coast where both models 
showed current wheat suitability, the GAM projected 
small yield losses (<10%; Fig. 5), further revealing the 
GAM'S tendency to show losses where DSSAT found 
gains. 

Our findings are consistent with EM-MM compari- 
sons conducted by ecological modelers. Morin & Thuil- 
ler (2009) used two EMs and two MMs to study North 
American tree species' range shifts in 2100, and found 
that the EMs projected larger range losses than the 
MMs. Cheaib et al.'s (2012) projections of Erench tree 
species distributions (for 2055) showed the same pattern. 
Although agricultural impacts researchers use both 
model classes, they are typically applied separately, and 
intercomparison studies focus on MM ensembles (Asseng 
et al., 2013; Rosenzweig et ai, 2013; Palosuo et ai, 2011; 
Challinor & Wheeler, 2008; Muller et al, 2011). However, 
separate studies suggest that EMs find larger yield losses 
than MMs (Schlenker & Lobell, 2010). 

There are several potential reasons why EMs tend to 
project larger climate losses than MMs. First, they 
assume that the modeled species currently occupies its 
full potential climate space, even though this may have 
been historically much broader (Morin & Thuiller, 


2009; Veloz et al., 2012). Regional agricultural impacts 
studies such as ours are sensitive to this problem. We 
relied on region-specific crop and climate data to fit 
models, even though these crops are grown globally 
under a broad range of climates. For example, if maize 
distribution data from neighboring countries were 
available to fit the GAM, it might have shown more tol- 
erance for the hotter, drier conditions along the margins 
(western through northeastern) of the current growing 
region (Figs 3 and 4). 

Second, EMs are generally fit to climatological vari- 
ables, but climate adaptation capacity may be better 
assessed by examining responses to interannual climate 
variability, rather than changes in climatic means. Eor 
instance, a tree species' climate tolerance may be best 
determined by the degree to which it can adapt phenol- 
ogy to interannual climate variations (Morin & Thuiller, 
2009). In agriculture, a crop's future viability may be 
more affected by failure rates related to changing inter- 
season variability rather than shifts in long-term means. 

Third, it is difficult for EMs to simulate management 
adaptations that can mitigate climate impacts. In our 
study, DSSAT models altered planting practices in 
response to both seasonal and long-term rainfall 
patterns, thereby simulating a certain level of manage- 
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men! adaptation that the GAMs could not. In this case, 
there was no apparent benefit, and the focus of our 
study was to understand between-model differences 
given current management and technology. This 
EM-MM difference would be more pronounced in a 
study explicitly focused on climate adaptation potential. 

Finally, as EMs do not simulate growth and repro- 
ductive processes, they cannot directly account for 
changes in the plant-climate relationship, such as 
increasing water-use efficiency and photosynthetic 
rates with rising CO 2 (Leakey, 2009). Cheaib et al. (2012) 
found CO 2 response to be a dominant factor influencing 
distribution projections for tree species whose range 
limits are determined by growth rate rather than cli- 
matic tolerances; projection discrepancies between EMs 
and MMs were particularly large for these species. 

Although EMs may overestimate climate-induced 
losses, MMs might underestimate them. One important 
reason relates to the aforementioned point on plant-C02 
interactions. The response of crops to elevated CO 2 is a 
key uncertainty (Parry et al., 2004; Tubiello et al., 2007), 
and some studies suggest that current MMs overesti- 
mate positive CO 2 effects (Ainsworth et al., 2008; Long 
et al., 2006). We found that elevated CO 2 increased 
DSSAT yield projections by 16% (maize) and 28% 
(wheat). DSSAT's CO 2 factor (Eqn 1) increases crop 
growth rates by 2% (maize) and 17% (wheat) for 
550 pm concentrations (relative to a 330 ppm baseline). 
At 660 ppm (100 + ppm greater than the levels we 
used) under well-watered conditions, DSSAT simulates 
ca. 4% maize yield increases (developer-reported values 
Hoogenboom et al., 2012), indicating that the CO 2 - 
related increases we found were primarily due to lower 
transpiration rates reducing water stress. There is no 
equivalent reported wheat CO 2 sensitivity, but our 
C02-related gain was 14% higher than the mean PACE- 
reported yield increase (Ainsworth et al., 2008), which 
suggests that this benefit was also due to increased 
water-use efficiency. We cannot judge how realistic 
these water-use efficiency gains are, but their size is a 
major source of difference between DSSAT and GAM 
projections, which had much smaller CO 2 adjustments 
applied. 

Crop MMs may also underestimate temperature 
sensitivity, particularly during flowering (Rotter et al., 
2011). Analyses suggest that recent warming has nega- 
tively impacted global maize and wheat yields (Lobell 
& Field, 2007), and EMs for major African crops 
(including maize) show much greater impacts from 
temperature than precipitation (Schlenker & Lobell, 
2010). African maize trial data show that temperatures 
>30 °C produce strongly negative and nonlinear 
impacts on yield (Lobell et al., 2011). In DSSAT, growth 
and grain filling halt at high temperatures, but sterility 


or mortality does not occur (Ritchie, 1998), thus temper- 
ature effects on crop reproduction may not be ade- 
quately represented. Low fertilizer parameter settings 
could also cause MMs to underestimate temperature 
sensitivity (Schlenker & Lobell, 2010). Although in this 
study, higher N inputs did not increase yield losses 
under warming (Appendix S4), the results also suggest 
that DSSAT was most sensitive to rainfall (e.g., very 
high predicted yields in high rainfall regions; Figs 2 
and 4), which was either constant or increasing under 
the climate scenarios we used. The GAMs may there- 
fore more accurately capture maize and wheat temper- 
ature sensitivities, as reflected by their yield loss 
projections even in areas of increasing rainfall (Figs 2, 
4, and 5). 

Mechanistic models could also underestimate nega- 
tive impacts by failing to simulate important biotic limi- 
tations (Morin & Thuiller, 2009), such as weed 
competition or pest damage. These are major crop 
growth impediments that are likely to be aggravated by 
climate change (Tubiello et al., 2007). DSSAT can simu- 
late pest damage, but we could not feasibly implement 
this for our study. In contrast, EMs implicitly account 
for these factors because they are encoded in the model 
training data (Morin &. Thuiller, 2009). The NDVIyieid 
proxies used to fit GAMs reflected not just climate 
and soil conditions but also pest prevalence, variable 
management practices, and other factors that determine 
actual productivity. Although the impacts of climate 
change on a crop's biotic interactions are likely to be 
complex, in some cases, competitors' or predators' 
niches may simply be shifted to new locations within 
the existing growing region, or have simple linear 
responses to climate change, in which case EMs may 
more realistically reflect potential impacts than MMs. 

Given the substantial differences between the two 
model classes' sensitivities, structures, and projected 
impacts, it is difficult to say that one class is inherently 
better than the other for studying climate change 
impacts. Instead, each may be better at simulating 
different aspects of crops' -climate responses (e.g., CO 2 
response by MMs, temperature sensitivity by EMs), 
while being sensitive to different input data errors, 
making it advantageous to compare EMs and MMs to 
obtain a fuller picture of projection uncertainties (Morin 
& Thuiller, 2009). 

Implications for South African cereal crops 

Our findings differ from previous studies providing 
nonspatial, country-level projections for South Africa, 
in that we find smaller negative impacts (GAM median 
—3.6%) or outright gains (DSSAT median -1-6.5%) for 
maize. In comparison, Schlenker & Lobell (2010) cited 
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Fig. 5 Decision Support System for Agrotechnology Transfer (DSSAT)- (top row) and Generalized Additive Model (GAM)-simulated 
(bottom row) wheat yields for the baseline (left column) and future (middle column) periods, and the differences (right column) 
between the two expressed as a percentage. 


three studies (including their own) projecting losses of 
10% or greater (a 4th showed yield increases, but was 
considered poorly calibrated). Fischer et al. (2005) in a 
study of cereal producfivity and suitability changes 
mapped projections for southern Africa showing severe 
productivity declines under 3 of 4 climafe scenarios 
within South Africa's maize production region. They 
also found substantial yield losses along the southern 
coast, where we either found small wheat yield losses 
(GAM) or substantial gains (DSSAT; Fig. 5). 

The difference between our results and the aforemen- 
tioned studies may be largely due to the CGCM-down- 
scaling technique we used (Hewitson &. Crane, 2006), 
which produced climate scenarios showing no change 
or modest precipitation gains over most of South Africa 
(Fig. 2), which contrasts with the losses projected by 
the IPCC AR4 for 2080 (IPCC, 2007). Although this 
downscaling method is well established (Flewitson & 
Crane, 2006; Patt et al., 2010; Tadross et al., 2009), 
understanding the food security implications of our 
findings in relation to previous estimates is difficult 
because earlier studies used substantially different 
climate scenarios. Further assessments are needed to 
test how different downscaling techniques (including 
both dynamic and empirical approaches) affect crop 
impacts projections for this region. 


Caveats 

Our study compared just one model per class (EM 
and MM) for each crop; therefore, the strength of our 
conclusions regarding systematic EM-MM differences 
is tempered by the small sample size. Flowever, this 
limitation is partially overcome by the consistency of 
our results from ecological model intercomparisons, 
and with the evident differences between EM and 
MM-based agricultural impacts studies. 

Soil data errors could also have biased our results. 
The soil data were of higher resolution than the global 
datasets used in comparable regional studies, but they 
only provided DSSAT's minimum parameter set, which 
were derived with simplified pedofransfer functions 
(Estes et al., 2013). Soil fertility was likely overesti- 
mated in high rainfall areas (Fig. 2), leading to inflated 
maize yields by DSSAT, particularly in the eastern 
highlands (Estes et al., 2013; Appendix S3), which is 
supported by the asymptotic relationship between 
NDVIyieid and rainfall in the GAM model (Appendix 
S3). DSSAT's seemingly greater sensitivity to rainfall 
than temperature may also have exaggerated this error. 
DSSAT also substantially underestimated yields in hea- 
vier clay soils due to overestimated WPs (Appendices 
SI and S3, Eig. 4). The combination of fhese fwo error 
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sources makes the DSSAT-NDVIyieid relationship 
appear to be nonlinear, but neither log transformation 
nor a second order polynomial model substantially 
improve the correlation between the two variables (Es- 
tes et al., 2013). Furthermore, the relationship between 
DSSAT and reported yields at the provincial scale is 
linear. 

The GAMs were less reliant on soil parameters and 
thus less prone to soil-related biases. The GAMs also 
were unaffected by the assumptions and errors inherent 
in choosing and calibrating the cultivar coefficients 
needed by DSSAT (Jones et al., 2003). We were unable 
to calibrate DSSAT for South African cultivars, which 
likely increased simulation error, as shown by a study 
comparing DSSAT and other MMs run with limited 
calibration (Palosuo et al., 2011). The greater sensitivity 
of MMs to data and calibration errors is another reason 
why it is valuable to compare results with EMs. 

The GAMs were affected by other weaknesses that 
could potentially alter interpretations. We followed the 
novel approach of fitting GAMs to a remotely derived 
yield proxy, rather than reported yields. A substantial 
body of work has found that NDVI-derived proxies accu- 
rately measure yields at a variety of spatial resolutions 
(see examples in Estes et al., 2013), which increases our 
confidence in this approach, but we could only validate 
this relationship at the provincial scale for maize. 

Another limitation was caused by the geographic 
distribution of the wheat dataset. As our GAM train- 
ing data were continuous observations (NDVIyiei^, as 
opposed to binary wheat presence/ absence data), the 
model could not be trained under climatically unsuit- 
able conditions because wheat is not grown in those 
areas. This limitation, combined with the small area of 
wheat cultivation in South Africa, made the model 
harder to constrain and less accurate than the maize 
GAM. Nevertheless, the wheat GAM was useful for 
comparing to DSSAT in their jointly suitable areas, 
where the relative differences in their yield projections 
reflected those seen between the maize models. 

Broader implications 

These findings suggest that climate impacts studies 
focused on food security and agro-ecosystems may be 
systematically influenced by the class of crop model 
used. Crop model intercomparison studies (e.g., 
Rosenzweig et al., 2013) should therefore expand 
beyond evaluating different MMs to include compari- 
sons between EMs and MMs, where differences may be 
larger. Because EMs and MMs characterize different 
aspects of the species-environment relationship (Mo- 
rin & Thuiller, 2009) and have different sensitivities, 
their combined use may provide a fuller picture of 


uncertainties in crops' climate responses than single- 
class ensembles. 
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